#ifndef HERMITE_INTERPOLATION_HPP  
#define HERMITE_INTERPOLATION_HPP  

#include <vector>
#include <functional>
#include <stdexcept>

class HermiteInterpolation {
public:
  
    // 默认构造函数
    HermiteInterpolation() = default;
    // 使用已知的数据点和导数值构造 Hermite 插值多项式
    HermiteInterpolation(const std::vector<double>& x_, const std::vector<double>& y_, const std::vector<double>& dy_) {
        if (x_.size() != y_.size() || x_.size() != dy_.size()) {
            throw std::invalid_argument("The sizes of x, y, and dy must be equal.");
        }
	
        int n = x_.size();
	x.resize(n);
	y.resize(n);
	dy.resize(n);

	n = 2*n;
	coef.resize(n);


        x.assign(x_.begin(), x_.end());
        y.assign(y_.begin(), y_.end());
	dy.assign(dy_.begin(), dy_.end());
	
	for (int i = 0; i < n; i++){
	  coef[i] = y[i/2];
	}
  
	for (int i = n-1; i >= 1; --i){
	  if (i%2 == 0)
	    coef[i] = coef[i] = (coef[i] - coef[i - 1]) / (x[i/2] - x[(i-1)/2]);
	  else
	    coef[i] = dy[i/2];
	}
  
	for (int level = 2; level < n; ++level) {
	  for (int i = n - 1; i >= level; --i) {
	    coef[i] = (coef[i] - coef[i - 1]) / (x[i/2] - x[(i-level)/2]);
	  }
	}
    }

    // 使用构造的 Hermite 插值多项式计算给定点的插值结果
    double interpolate(double z) const {
	int n = coef.size();
        double result = coef[n-1];
        for (int i = static_cast<int>(coef.size()) - 2; i >= 0; i--) {
            result = result * (z - x[i/2]) + coef[i];
        }
        return result;
  };
    

private:
    std::vector<double> coef; // 存储插值多项式的系数
    std::vector<double> x;    // 存储原始的 x 值，用于插值时计算
    std::vector<double> y;    // 存储原始的 y 值
    std::vector<double> dy;   // 存储原始的导数值
};

#endif
